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Detailed protocol 


Dear Yankai, 

Thank you for the nice words on our paper. From what | understood your have question regarding sample integration. We actually used the Seurat v3 pipeline 
which relies on anchor cells to compute batch-correction vectors between datasets. https://satijalab.org/seurat/archive/v3.0/integration.html| 

You can find the Seurat code below. | have also sent you this script via email since bioprotocol doesn't allow to send R scripts. 

Best, 


library(Seurat) 
library(gplots) 
library(paletteer) 
library(tidyverse) 
# First we generate one seurat object for each patient and put them ina list 
eleven_cd3.list <- c(P34_Tumor,P35_Tumor,P42_Tumor,P43_Tumor,P46_Tumor,P47_Tumor,P55_ Tumor, P57_Tumor, P58_ Tumor, P60_Tumor, 
P61_ Tumor) 
for (i in 1:length(x = eleven_cd3.list)) { 
eleven_cd3 list{[i]] <- NormalizeData(object = eleven_cd3.list[[i]], verbose = T) 
eleven_cd3 list{[i]] <- FindVariableFeatures(object = eleven_cd3.list[[i]], 
selection.method = "vst", nfeatures = 2000, verbose = T) 

} 
eleven_cd3.anchors <- FindIntegrationAnchors(object.list = eleven_cd3.list, dims = 1:30) 
eleven.tils.cd3.integrated <- IntegrateData(anchorset = eleven_cd3.anchors, dims = 1:30) 
# switch to integrated assay. The variable features of this assay are automatically set during IntegrateData 
DefaultAssay(object = eleven.tils.cd3.integrated) <- "integrated" 
# Run the standard workflow for visualization and clustering 
eleven.tils.cd3.integrated <- ScaleData(object = eleven.tils.cd3.integrated, verbose = FALSE) 
eleven.tils.cd3.integrated <- RunPCA(object = eleven.tils.cd3.integrated, npcs = 30, verbose = FALSE) 
eleven.tils.cd3.integrated <- FindNeighbors(object = eleven.tils.cd3.integrated, dims = 1:30) 
eleven.tils.cd3.integrated <- FindClusters(object = eleven.tils.cd3.integrated, resolution = 1.2) 
eleven.tils.cd3.integrated <- RunUMAP (object = eleven.tils.cd3.integrated, reduction = "pca", 

dims = 1:30) 
DimPlot(object = eleven.tils.cd3.integrated, reduction = "umap", group.by = "orig.ident", pt.size = 1) 
DimPlot(object = eleven.tils.cd3.integrated, reduction = "umap", pt.size = 1, label = T, repel = T) 
FeaturePlot(eleven.tils.cd3.integrated, "CD3D", pt.size = 1, min.cutoff = "q3" order = T, max.cutoff = "q97") 
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